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Abstract: The shrinkage of the Aral Sea, which is closely related to the Amu Darya River, strongly affects 
the sustainability of the local natural ecosystem, agricultural production, and human well-being. In this 
study, we used the Bayesian Estimator of Abrupt change, Seasonal change, and Trend (BEAST) model to 
detect the historical change points in the variation of the Aral Sea and the Amu Darya River and analyse 
the causes of the Aral Sea shrinkage during the 1950-2016 period. Further, we applied multifractal 
detrend cross-correlation analysis (MF-DCCA) and quantitative analysis to investigate the responses of 
the Aral Sea to the runoff in the Amu Darya River, which is the main source of recharge to the Aral Sea. 
Our results showed that two significant trend change points in the water volume change of the Aral Sea 
occurred, in 1961 and 1974. Before 1961, the water volume in the Aral Sea was stable, after which it began 
to shrink, with a shrinkage rate fluctuating around 15.21 km/a. After 1974, the water volume of the Aral 
Sea decreased substantially at a rate of up to 48.97 km3/a, which was the highest value recorded in this 
study. In addition, although the response of the Aral Sea's water volume to its recharge runoff 
demonstrated a complex non-linear relationship, the replenishment of the Aral Sea by the runoff in the 
lower reaches of the Amu Darya River was identified as the dominant factor affecting the Aral Sea 
shrinkage. Based on the scenario analyses, we concluded that it is possible to slow down the retreat of the 
Aral Sea and restore its ecosystem by increasing the efficiency of agricultural water use, decreasing 
agricultural water use in the middle and lower reaches, reducing ineffective evaporation from reservoirs 
and wetlands, and increasing the water coming from the lower reaches of the Amu Darya River to the 
1961-1973 level. These measures would maintain and stabilise the water area and water volume of the 
Aral Sea in a state of ecological restoration. Therefore, this study focuses on how human consumption of 
recharge runoff affects the Aral Sea and provides scientific perspective on its ecological conservation and 
sustainable development. 
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1 Introduction 


Water is an important environmental factor and natural resource that plays a crucial role in natural 
ecosystems, agricultural production, and human well-being. However, freshwater systems, 
including rivers, lakes, and springs, are directly or indirectly subjected to human influence 
through land cover change and anthropogenic climate change (Chen et al., 2019; Konapala et al., 
2020). This has led to severe consequences, such as the alteration of waterways and the 
desiccation of rivers and lakes, resulting in ecological and environmental problems including 
water pollution, ecosystem destruction, habitat loss, etc. For instance, studies have reported that 
approximately 20% of freshwater lakes have become desiccated, while the altered volume has 
decreased by over 40% during the last century (Feng et al., 2022), causing 18% of the world's 
population to face water shortage challenges (Berdugo et al., 2022). Considering the rapid 
expansion of global drylands (Huang et al., 2016) and population growth (Frona et al., 2021), 
these problems are likely to continue and even be exacerbated in the future (Wang et al., 2018; 
Woolway et al., 2020; Grant L et al., 2021). Because of their severity, the sustainable use of water 
resources has been added to the UN sustainable development goals (Foster, 2018). As water 
resources are vulnerable to the harmful impacts of climate change, water shortage and desiccation 
are especially more pronounced in arid lands owing to low per capita natural resources and 
inadequate basic infrastructure and water-saving technology in agriculture. 

Central Asian is the world's largest arid region, of which the water system mainly includes 
glacier-fed surface runoff and groundwater. Owing to the intensive overexploitation of water 
resources and the increase in evaporation and rapid glacier retreat during the past century due to 
climate change, Central Asia has suffered from severe water shortage, leading to serious 
environmental problems, including desert expansion, soil salinization, sand storm occurrence, and 
water pollution (Zhang et al., 2019; Pham-Duc et al., 2020; Hu et al., 2022). Consequently, 
according to a recent report by Li et al. (2022), soil salinization and desertification-related 
economic losses are as high as 3.0x10'° USD, posing a serious threat to the food security and 
sustainable development of the region. 

Inland water systems, such as terminal lakes, are essential for maintaining ecosystem stability 
and regional safety. The desiccation of lakes has a detrimental effect on the natural environment 
and local population. For example, the Lop Nur Lake desiccation has led to a substantial 
expansion of the sandy desert, the disappearance of small oases, the complete destruction of the 
regional ecosystem, and severe soil salinization (Cai et al., 2022). To prevent similar catastrophes, 
understanding lake distribution and variation in Central Asia and their underlying mechanisms has 
been a common area of interest in climatic, hydrologic, ecologic, and geographic research (Bai et 
al., 2012; Liu et al., 2019). As the largest inland terminal lake in arid Central Asia, the Aral Sea 
plays a key role in the ecological environment and water resource utilisation. However, along 
with population growth, the rapid expansion of croplands in the middle and upper reaches of the 
Amu Darya River and Syr Darya River (Rufin et al., 2022), which represent the main water 
sources of the Aral Sea, has significantly decreased the incoming water to the lake basin. In 
addition, over the past 100 years, glacial recession and an increase in evapotranspiration caused 
by global warming have exacerbated the severity of this problem in Central Asia (Wang et al., 
2020). Furthermore, the construction of hydrological infrastructure, such as hydropower stations 
built in Tajikistan, Uzbekistan, and Kyrgyzstan in the upper reaches of the Amu Darya River and 
Syr Darya River, and water drainage channels, including the Karshi Canal, Karakum Channel, and 
Amu-Bukhara Canal, have also reduced the water flow into the lake. Moreover, a lack of unified 
systematic monitoring, management, and planning measures for water resources and agricultural 
lands—especially after the collapse of the Soviet Union—has led to excessive land reclamation, 
overexploitation of groundwater, and inappropriate irrigation, causing a reduction in the runoff 
volume of rivers, the decline of the water table, and secondary salinization (Liu et al., 2022). 
Consequently, in 1951, 6.600x10* km? of the water surface area of the Aral Sea, which was 
previously used to support local fishery, was divided into four separate water patches spanning 
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7.000x10° km? today (Wu et al., 2022), creating large swathes of saline islands and sandy deserts 
around and in the middle of the lake basin. This has led to the desiccation of the local climate, 
causing frequent heatwaves and sandstorms that pose serious threats to socioeconomic 
development and even the health and well-being of the local population (Sharma et al., 2018; He 
et al., 2022). 

Thus, understanding the hydrological and anthropogenic processes that caused lake desiccation 
in order to find a way to mitigate it is one of the main focuses of dryland research. For instance, 
research conducted by Huang et al. (2022) and Zhang et al. (2022), who analysed glacier changes 
over the past century, showed that the inland lakes in Central Asia decreased, while the 
higher-altitude alpine lakes increased significantly along with glacier meltwater, stressing the 
importance of human activities in the lake area reduction of the Aral Sea. Rufin et al. (2022) 
reported that croplands increased from 4.400x10* km? in 1987 to 5.200x10* km? in 2017, with 
more than a 200% increase in double-cropping irrigation lands, followed by a rapid urban 
expansion and sharp rise in water consumption in the last 20 years (Yan and Li, 2023), 
contributing to the problem. A comparative analysis of climatic and anthropogenic factors by 
Yang et al. (2020a) concluded that anthropogenic factors were the main cause of the reduction in 
the lake area. Even though these studies have deepened our understanding of the relative 
contributions of climatic and anthropogenic factors to some extent, the ways in which specific 
hydrological and human-driven processes of the major water sources of the lake have led to the 
lake area reduction of the Aral Sea remain unclear. As a major water source, the Amu Darya River 
has considerable influence on the lake area fluctuations. Human activities or water volume 
fluctuations in any reach of the Amu Darya River could directly affect the lake area variation. 
However, the effect of runoff variations in the Amu Darya River on Aral Sea water area 
fluctuations is still mostly unknown. Therefore, analysing the response of the Aral Sea's water 
volume and lake area to the runoff variations of the Amu Darya River would help us to 
understand the mechanisms behind the shrinkage of the Aral Sea. Further, it would provide a 
scientific basis for local authorities to improve their water resource management regime and halt 
the continued devastation of the region's natural environment. 

In general, parametric (e.g., cumulative anomaly and sliding t-test) and non-parametric (e.g., 
Mann-Kendall test, Pettitt's test, and Bayesian analysis) methods (Hrvatin and Zorn, 2022) are 
widely used for detecting the change-point in surface runoff as well as the corresponding water 
volume and lake area. These methods are easy to operate and can deliver relatively robust results 
when observation data are sufficient and the relationship between the explanatory factors of the 
lake area and volume variation appears to be relatively simple (Alkan and Konukcu, 2022). As 
most of these methods are based on statistical parameters, such as Akaike's information criterion 
and the Bayesian information criterion (Hrvatin and Zorn, 2022), their produced results do not 
always align with each other and can even be contradictory (Zhao et al., 2019). Conventional 
methods may not necessarily perform well when analysing complex hydrological processes. A 
holistic approach based on the direct and indirect hydrological processes of river networks has 
been proven to more realistically reflect the underlying mechanisms of lake variation (Zhao et al., 
2004). In this study, we used a novel approach for hydrological analysis—the Bayesian Estimator 
of Abrupt change, Seasonal change, and Trend (BEAST) model—to explore this relationship. As 
the BEAST model comprehensively ensembles various competing models based on the Bayesian 
approach, it can effectively help explore intrinsic causal relations in terms of the dimension of 
temporal dynamics (Zhao et al., 2019). Other popular change-point detection methods are based 
on other single values, such as the mean and variance (Pettitt, 1979; Vinushree et al., 2022; 
Banakara et al., 2023), which are often unable to precisely find change-points, leading to a 
contradiction in the explanatory factors. 

Furthermore, we used the multifractal detrend cross-correlation analysis (MF-DCCA) method 
to explore the complex non-linear relationships between the hydrological variables in this study; 
the MF-DCCA method shows some advantages in handling the non-linear interaction relations 
between non-stationary time series, which has been widely applied in finance but rarely used in 
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hydrology. This method discerns the non-linear dynamic pattern of each variable by both studying 
the multifractal features of the variable sequence and calculating the long-range correlation of the 
variable obtained by varying the time-scale of the Hurst index (Xiang et al., 2019; Rahmani and 
Fattahi, 2021; Yang et al., 2021; Kojić et al., 2022). 

Based on these considerations, the surface runoff data from the hydrological monitoring 
stations along the Amu Darya River, and the satellite-driven water level, water area, and water 
volume data of the Aral Sea from 1950 to 2016, we analysed the effect of Amu Darya River 
runoff variations on the Aral Sea fluctuations to provide new insights and a scientific basis for 
sustainable water use and management in the region. Specifically, we addressed the following 
questions: 

(1) What are the dynamics of the Aral Sea and the surface runoff and water consumption of the 
Amu Darya River during 1950-2016? 

(2) What are the reasons for the Aral Sea shrinkage, and do they vary with time? 

(3) How does the Amu Darya River runoff affect the water volume of the Aral Sea during the 
study period? 

Analysing the connection between the Aral Sea shrinkage and the Amu Darya River will not 
only clarify the role of human activities in the evolution of rivers and lakes in the arid zone of 
Central Asia, but also provide a scientific basis for the rational development and protection of 
water resources and the promotion of regional sustainable development. 


2 Data and methods 


2.1 Study area 


The Aral Sea (58°-62°E, 43°-47°N) is located in Central Asia, at the border of the Republic of 
Karakalpakstan in Uzbekistan and the Kyzylorda and Aktobe oblasts in Kazakhstan. It meets the 
delta of the Syr Darya River to the north and the delta of the Amu Darya River to the south, with a 
high topography in the west and a low one in the east (Fig. 1). The Aral Sea region is 
characterised by an extreme continental climate. Precipitation is higher in April and October (with 
a daily amount of 0.48 mm) and lower in August and September (with a daily amount of 
approximately 0.22 mm). Temperature is higher in summer (with a monthly average value of 
25°C and a monthly maximum value of 47°C) and lower in winter (with a monthly average 
temperature of —11°C and a monthly minimum temperature of —38°C). Annual evaporation is as 
high as 1700 mm. As a terminal lake in an arid zone in Central Asia, runoff from the Amu Darya 
River and Syr Darya River is the largest source of its recharge. 

The Aral Sea was once the second largest inland saltwater lake in Asia after the Caspian Sea 
and the fourth largest lake in the world. In the 1950s, the Aral Sea had a water area over 
7.000*10* km, an average altitude of 53 m, a maximum length of 435 km from north to south, a 
width of 290 km from east to west, an average depth of 16 m, a maximum depth of 69 m on the 
west bank, and a total water volume of approximately 1.000x10* km? (Micklin, 1988, 2007, 2010, 
2016; Micklin and Aladin, 2008). In the 1960s, the Aral Sea covered an area of approximately 
6.500x10* km’, but its area decreased year by year, splitting into the South Aral Sea (Big Aral 
Sea) and North Aral Sea (Small Aral Sea) in 1986, with the Big Aral Sea further splitting into the 
East Aral Sea and West Aral Sea in 2007. The lake area has now been reduced to approximately 
0.400x10* km’, with the elevation of the lake surface in the West Aral Sea dropping to 26 m, 
while the Small Aral Sea remains above 40 m, and the total water volume of the lake has dropped 
to about 50.000 km. Overall, the Aral Sea has shrunk rapidly since the 1960s, with its lake area 
now being less than 10% of its previous size, leading to huge impacts on the local ecological 
environment and social economy. 


2.2 Data sources 


The data used in this study mainly include lake data from the Aral Sea (before 1986) and Big Aral 
Sea (after 1986), hydrological data from the Amu Darya River and Syr Darya River, and 
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Fig. 1 Overview of the Aral Sea region in Central Asia (a) and MODIS images showing the lake area variations 
of the Aral Sea in 1977 (b), 2010 (c), 2020 (d). MODIS images are derived from National Aeronautics and Space 
Administration (NASA) (https://www.earthdata.nasa.gov/sensors/modis). 


socioeconomic data from the Aral Sea, which are all annual data. The lake data, obtained from the 
CA Water website (http://www.cawater-info.net/), comprise the water level, water area, and water 
volume during the 1950-2017 period and the theoretical capacity of the Aral Sea. 

The hydrological data include the in-situ runoff data (1950-2016) from four hydrological 
stations located in the upper, middle, and lower reaches of the Amu Darya River (Kelif, Kerki, 
and Samanbay hydrological stations, respectively) and the lower reaches of the Syr Darya River 
(Karateren hydrological station, denoted by the symbol Syr_ Down). Additionally, they include the 
reservoir storage capacity (1983-2017) of the Tuyamuyun Hydro Complex (THC), which is 
located at the beginning of the lower Amu Darya River, between the Kerki and Samanbay 
hydrological stations, and is the sole and largest reservoir in the downstream area, with a total 
capacity of 6.900-7.800 km?. 


2.3 Correction of the lake data 
Because of inconsistencies in the original historical annual data, their direct use leads to distorted 
results and unsatisfactory simulations. These inconsistencies are mainly rooted in two aspects: (1) 


the correspondence between the three types of lake data (water level, water area, and water 
volume) fails to follow the theoretical capacity curve and (2) the inconsistency between these lake 
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datasets in their temporal fluctuations. This is attributed to the inconsistency of temporal data we 
obtained, as the Aral Sea fluctuated greatly within an annual scale, and its water volume values 
varied greatly among different months. In addition, the raw data were prone to errors in the 
process of being recorded, and the raw data were not collected at the same or similar times. 
Therefore, we developed an analytical procedure to derive the most suitable time series data using 
a combination of mapping based on theoretical capacity curves and optimisation techniques, as 
described in Figure 2. 
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Fig. 2 Flow chart of the lake data correction in this study. WV, WL, and WA represent the water volume, water 
level, and water area, respectively; WV_2 and WV_3 represent the transformed water volume under the two 
transformation methods, respectively; WA_c, WL_c, and WVC ¢ represent the corrected water area, water level, 
and water volume change, respectively; H-W process and S—W process represent the inversion of the water 
level—water volume relation and the inversion of the water area—water volume relation, respectively. 


In this study, we focused on the water volume data and its changes. The sets of original annual 
water level and water area data for the Big Aral Sea from 1950 to 2017 were first mapped into the 
water volume data using the theoretical capacity curves. Through mapping, we obtained three 
types of the water volume data: the original data, the transformed data obtained by the inversion 
of the water level—water volume relation (H—W process), and the transformed data obtained by 
the inversion of the water area—water volume relation (S—W process). The Aral Sea was split into 
the Small Aral Sea and Big Aral Sea in 1986, while the latter was further split into the East Aral 
Sea and West Aral Sea in 2007. We fit the water level-water area—water volume relationship with 
polynomials and considered those in three segments corresponding to three periods in history: 
1950-1985, 1986-2006, and 2007-2017. Finally, we obtained several sets of water volume data 
from 1950 to 2017, constituting a range of values. We selected the most suitable single values for 
the follow-up study from the range of values constituted by the available water volume data using 
optimisation techniques, whose essence lies in the temporal inconsistency of the data. As it is 
difficult to obtain long time series data that are mostly consistent over time, optimisation 
techniques are useful and necessary. During the optimisation, the constrained non-linear 
minimisation solver was used to minimise the negative correlation between the water volume 
change and its recharge runoff. Using the interior point algorithm and 355 iterations, we obtained 
the final water volume data. As for the runoff recharge, the downstream recharge of the Amu 
Darya River and Syr Darya River was combined before 1986, after which only the downstream 
recharge of the Amu Darya River was considered, with the corresponding water volume change 
from the Big Aral Sea. As the shrinkage of the Aral Sea mainly occurred in the Big Aral Sea 
during 1986-2017, the Small Aral Sea was isolated from the Big Aral Sea, with the water volume 
increasing steadily. Using the theoretical capacity curves, the corresponding corrected water level 
and water area data were obtained by the inversion of the corrected water volume data. 


2.4 Sen's slope calculation 


Sen's slope calculates the slope (i.e., the linear rate of change) and confidence level (Pohlert, 
2016), according to the following equations: 
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where dj is the slope calculated between the two time points t; and 4; x; and x; are the 
corresponding values of the variables at time points t; and tj, respectively; N is the length of the 
data; bsen is the Sen's slope; and d is the set of all slopes with a total number of M(N-1)/2. 


2.5 Multifractal detrend cross-correlation analysis (MF-DCCA) method 


The MF-DCCA method proposed by Zhou (2008) is based on the multifractal detrended 
fluctuation analysis and detrended cross-correlation analysis (DCCA) methods. It mainly includes 
the following steps. First, by assuming that two non-stationary sequences x(u) and y(u) have the 
same length N (where u=1, 2, ..., N), we constructed two new sequences X(/) and Y(/) (1, 2, ..., 
N) according to Equation 3: 


X(N=DVi(x()-+')s YO= VOW) -y'), 6) 
u=l u=l 

where / is the index for each variable in the new sequences X and Y; and x’ and y’ represent the 
mean values of the time series x(u) and y(u), respectively. 

Next, the new sequences X(/) and Y(/) were divided into N, equal-length and non-overlapping 
sub-intervals, with Ns=int[N/s], where s represents the length of the sub-intervals (i.e., timescale) 
and N; is the total number of sub-intervals. Because the timescale s may not be an integer multiple 
of N, the time series were again segmented from back to forward, with 2; sub-intervals for every 
time series. To obtain the fluctuation function of order g, we averaged the detrended variance for 
all sub-intervals. When g=0, the formula was derived from L'H6pital's rule: 


1/q 
F,(s)= lade [F (s, ia 4+0, () 
F, ()=e| ye In| F? (=»)]} q=0, (5) 


where F(s) is the fluctuation function of order q when the sub-interval's length is s; Fo(s) is the 
fluctuation function when g=0; v is the index for each sub-interval; and F?(s,v) is the detrended 
variance function of each sub-interval, which is computed by calculating the variance removing 
the local trend obtained by the least-squares fitting of each sub-interval, as shown in Equations 4 
and 5. When v=1, 2, 3, ..., Ns, then 


F*(s,v)=2 SA x[(v-s+e]-%, (eV [ols +2 ]-¥ (2) 


S 


2 (6) 


g= 


and when v=N,+1, N.+2, Ns+3, ..., 2s, then 


Psa DAN- 0M) +] (PLY (v-N,)s+z]-¥, (r)}. (7) 


Here, X\(z) and Y,(z) are the local trends in each sub-interval, and z is the index of each variable in 
the sub-interval (z =1, 2, ..., s). We repeated the described steps using different values of g and 
examined how F;,(s) changed with timescale s by observing the double log plot between F;,(s) and 
s. When a power-law cross-correlation exists, the relation satisfies the following formula: 


F (s) ss sho) : (8) 


q 
where s**»®) represents the function of the timescale s, and Ax,(q) is a specific exponent called the 
generalized Hurst exponent. When g=2, the MF-DCCA becomes the DCCA. If hx,(2)=0.5, no 
cross-correlation exists between the two time series. If h,,(2)>0.5, then a forward long-range 
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cross-correlation exists; whereas if /,,(2)<0.5, an inverted long-range cross-correlation is 
observed between the two time series. The multifractal features of a time series can be expressed 
by the dependence of x(q) on q, whereas the monofractal data h,,(q) are independent of the 
relative q values. When g>0, the generalised Hurst exponent h,,(q) describes the scale effect of 
large fluctuations, whereas when q<0, /x;(qg) characterizes the scale effect of small fluctuations. 


2.6 Bayesian Estimator of Abrupt change, Seasonal change, and Trend (BEAST) model 


The BEAST model was developed to isolate periodic and trend signals from a time series and 
pinpoint abrupt shifts by fitting a Bayesian regression model to match the time series data (Zhao 
et al., 2019). A time series D={tu, Yu} (where t, is the time of each observation; y, is the observed 
value at time ¢,; and u is the index for each observation (u=1, 2, ..., N)) is believed to comprise 
three components (seasonality, trend, and mutation plus noise), which can be expressed using the 
following formula: 

Vy =S(t,;0;)+T(t,;0,)+6&,. (9) 
Here, the noise €, is assumed to be Gaussian, with a magnitude o, and the rest of the data are 
explained by the seasonal S(¢,; Qs) and trend 7(t,; Or) signals (where Os is the parameter of the 
seasonal signals; and Or is the parameter of the trend signals). The general linear model was 
utilised to parameterise S(t,; Os) and 7(t,; Or). The mutation point (i.e., change point) is implied 
in the Os and Or parameters. 

For formulation purposes, the Or and Os parameters were reclassified into two sets, i.e., {O7, 
Os\}={M, Bu}. The first set, M, refers to the model structure, including the number and time of the 
seasonal and trend change points and the seasonal harmonic order. The second set, fy, is the 
coefficient parameter of a particular segment used to determine the exact trend shape and seasonal 
curves after a given model structure M. After recombination, the original general linear model 
formula becomes: 

y(t,) =Xu (ti Bu +> (10) 
where y(t.) is the observed value at time ¢,; and Xi(t,) is the dependent variable at time tu. This 
equation was extended to a general linear model to create a Bayesian model used to detect sudden 
changes, seasonality, and trends from the time series data. In the Bayesian modelling, all 
unknown parameters were considered as random, including the model structure M, the coefficient 
Bu, and the data noise o°. In a given time series D={t,, yu}, it is necessary to both calculate the 
best parameter values and obtain their posterior probability distribution p(B, o°, M|D), which is 
the product of a likelihood function and a prior model: 


P(By.0"M |D) x p(D| Buo”, M)xa( Buo, M), (11) 


where D is the observed data; p(D|Bm, o’, M) is the probability of observing the data D given the 
model parameters fu, o’, and M, and z(fu, o’, M) is the prior distribution of the model 
parameters. 

The posterior probability distribution p(B, o°, M|D) contains all the necessary information to 
infer the dynamic features of the system, including trends, seasonal changes, and mutations. The 
Markov chain Monte Carlo (MCMC) sampling method was used to achieve posterior inference. 
The study utilised a mixed sampling algorithm to embed the reverse-jump MCMC sampler into 
the Gibbs' sampling framework, as detailed in Zhao et al. (2019). 


2.7 Quantitative analysis of the impact of the Amu Darya River on the Aral Sea's water 
volume 


To quantify the runoff effects of the Amu Darya River on the water volume of the Big Aral Sea, 
we built a model according to a combination of the water balance equation, multiple regression, 
and polynomial fitting. First, based on the water balance equation and the closed characteristics of 
the Aral Sea region, we assumed that the amount of the water volume change equalled to the total 
evapotranspiration minus the total recharge, as expressed in Equation 12: 


WVC=(Eva — Pre) x WA — g(R) , (12) 
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where WVC is the water volume change (km?/a); Eva is the evaporation (mm); Pre is the 
precipitation (mm); and WA represents the water area (km7), and when WA is multiplied by the 
difference between evaporation and precipitation, the net evaporation (i.e., the expense item) can 
be obtained. The income term, including various forms of runoff recharge distributed in different 
spatial locations, can then be considered a function of the runoff of the supply river, defined by 
g(R). Water volume change, which is equal to the difference between the water volume at a given 
time m (WV »; km?) and the water volume at a subsequent time m+1 (WV m+1; km°), represents the 
water volume variation of the lake body; the (—) sign indicates a decreasing water volume (i.e., 
shrinkage), while the (+) sign indicates an increasing water volume (i.e., expansion). 

When the difference between evaporation and precipitation is assumed to be constant or have a 
small variation and the efficiency of the runoff recharge is constant over a certain time range, the 
calculation of the magnitude of the water volume change can be converted into a linear simulation 
of variables water area and recharge runoff from river n (Rn; km*) using the multiple regression 
method. Equation 12 was, thus, converted into Equation 13, and the corresponding coefficients a 
and p could be estimated using the least squares method: 


WVC=axWA-)) 8, xR, . (13) 
Furthermore, water area is a function of water volume, also called storage, as expressed in 
Equation 14. Because they are determined by the topographic features of the same lake basin, 


there exists a one-to-one correspondence between the water area and water volume. To continue 
combining Equations 13 and 14, we obtained the piecewise function in Equation 15: 


a,x f,(WV)- BR, — B,xR,, WV e WV, 
a, x f,(WV)-B, xR, — B,xR,, WV e WV, 
a,x f,(WV)-By, xR, WV e WV, 
a,x f,(WV)- 8, xR, WV e WV, 


WVC= (15) 


where Rı represents the recharge runoff from the Amu Darya River (km°), for which runoff data 
from the Samanbay hydrological station were used, while R2 denotes the recharge runoff from the 
Syr Darya River (km°), for which the runoff data from the Karateren hydrological station were 
utilised. There are four main segments in Equation 15 that are determined by three important time 
points (1966, 1986, and 2007) in the historical period when the Aral Sea underwent dramatic 
changes. Before 1986, the Aral Sea was simultaneously recharged by the Amu Darya River and 
Syr Darya River; however, after 1986, the Aral Sea receded into the Big Aral Sea, which was 
recharged almost exclusively by the Amu Darya River. Thus, WV;, WV2, WV3, and WV4 
represent the water volume separately during the periods of 1950-1965, 1966-1985, 1986-2006, 
and 2007-2016; fi, 4, f, and f4 represent the function describing the correspondence between the 
water area and water volume separately during the four periods; and a1, a2, Bi, £12, and f2 are the 
coefficients. 

The downstream runoff from the Amu Darya River can also be reflected by the water used in 
the middle-upper reaches (Cı; km*), water used in the middle-lower reaches (C2; km°), 
upstream-originating flow (U; km*), and water storage in the reservoir of the THC (RV; km°), as 
described in Equation 16. To continue studying the influence of factors, such as the water used in 
each river section of the Amu Darya River, on the Big Aral Sea, we constructed piecewise 
Equation 17. Owing to the properties of the multiple regression method, we chose the used water 
as independent variable instead of runoff whenever possible, to avoid multi-collinearity. 

R =h(C,,C,,U,RV), (16) 
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œ x fh (WV) -4 XG -7 XC =7; xU -y4 x RV, WV € WV; 

oœ x h (WV)-7 xC -7 xC -7 xU — y, x RV, WV e WV,” 
where a3, yı, 2, y3, and y4 are the coefficients. Here, the formula is divided into two segments 
corresponding to the stages before and after the Big Aral Sea split in 2007. As the THC reservoir 


was built in 1983, the water storage in the reservoir of the THC operated throughout the Big Aral 
Sea period. 


2.8 Set of variables 


Characterising the temporal dynamics of the water volume change may not only help understand 
the manner in which the water volume in the Aral Sea changed during the study period but also 
help explore the relationship between the water volume and its recharge runoff more directly and 
effectively in terms of the water balance equation. Therefore, this study chose the water volume 
change and the following variables related to its recharge runoff for this study: the runoff from the 
downstream of the Syr Darya River (Syr_Down_R), the runoff from the downstream of the Amu 
Darya River (Samanbay_R), the water consumption in the middle-upper reaches of the Amu 
Darya River (obtained by subtracting the runoff at the Kelif hydrological station from Kerki 
hydrological station (K-K)), the water consumption in the middle-lower reaches of the Amu 
Darya River (obtained by subtracting the runoff at the Kerki hydrological station from Samanbay 
hydrological station (K-S)), the water storage in the THC (THC_V), the runoff from the 
midstream of the Amu Darya River (Kerki_R), and the upstream-originating flow of the Amu 
Darya River (Kelif_R). Compared to other variables located in the Amu Darya River (from 1950 
to 2016), the time series of Syr Down_R was shorter (from 1950 to 1985) owing to the historical 
context of its recharge into the Big Aral Sea, which is also explained in the following quantitative 
analysis. 


wye=| (17) 


3 Results 


3.1 Dynamics of the water volume change in the Aral Sea and its recharge runoff 


3.1.1 Change point detection based on the BEAST model 


The trend change point detection results for each variable (water volume change, Samanbay_R, 
K-S, K-K, Kelif R, and Syr Down_R) provided the number and year of change points as 
probabilities. The model results showed the probability of the number of change points for 
different variables. First, we selected all the numbers that had a probability greater than 0.100 
(P>0.100). The water volume change showed the highest probability when it had 2 and 3 change 
points, with values equalling 0.532 and 0.326, respectively, and the 3 change points appeared in 
1961 (P=0.725), 1974 (P=0.932), and 1992 (P=0.142). The Samanbay_R variable showed the 
highest probability when it had 1, 2, and 3 change points, with values equalling 0.387, 0.380, and 
0.180, respectively, and the 3 change points occurred in 1961 (P=0.298), 1974 (P=0.939) and 
1969 (P=0.074). The K-S variable also showed the highest probability when it had 1, 2, and 3 
change points, with values of 0.397, 0.376, and 0.177, respectively, and the 3 change points 
appeared in 1961 (P=0.284), 1974 (P=0.931), and 1969 (P=0.073). The K-K variable had its 
highest probability with 2 and 3 change points, showing values of 0.856 and 0.101, respectively, 
and the 3 change points appeared in 1959 (P=0.953), 1981 (P=0.175), and 1991 (P=0.197). The 
Kelif_R variable exhibited the highest probability when it had 0, 1, and 2 change points, with 
values equalling 0.345, 0.370, and 0.159, respectively, and the 2 change points occurred in 1969 
(P=0.084) and 1974 (P=0.106). Finally, the Syr_ Down_R variable had the highest probability at 
3 and 4 change points, with values of 0.774 and 0.210, respectively, and the 4 change points took 
place in 1957 (P=0.819), 1961 (P=0.991), 1965 (P=0.154), and 1974 (P=0.802). We further 
observed the change point years of each variable obtained from the initial screening, then 
eliminated the years with significantly lower probabilities, selected the years with significantly 
higher probabilities, and plotted them in Figure 3 (dashed line). 
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Fig. 3 Results of the change point detection for the study variables based on the Bayesian Estimator of Abrupt 
change, Seasonal change, and Trend (BEAST) model. (a), water volume change; (b), Samanbay_R (runoff from 
the downstream of the Amu Darya River); (c), K-S, water consumption in the middle-lower reaches of the Amu 
Darya River; (d), K-K, water consumption in the middle-upper reaches of the Amu Darya River; (e), Kelif_R, 
upstream-originating flow of the Amu Darya River; (f), Syr_Down_R, runoff from the downstream of the Syr 
Darya River. Dashed lines represent the location of change points, and grey envelopes indicate 95% credible 
intervals for the fitted trend signals. 


The Samanbay_R and Syr_Down_R variables, representing the recharge runoff of natural 
channels in the Aral Sea, for the Amu Darya River and Syr Darya River, respectively, had the 
same change points as the water volume change in 1961 and 1974, showing a strong connection. 
The K-S and Samanbay_R variables showed high consistency in the change point detection 
results, denoting the high intensity of the connection. However, the results of the K-K variable 
completely differed from those of the Samanbay_R variable, and also did not present the same 
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change point years as the water volume change. This indicates that, compared with the K-S 
variable, the K-K variable had a smaller effect on the downstream runoff of the Amu Darya River 
and a weaker direct connection to the Aral Sea, suggesting that there may be more complex 
interactions; further argumentation may be given in the future. The Kelif_R variable had a change 
point year in 1974, as did the Samanbay_R variable and water volume change, but its occurrence 
probability value was rather low (only 0.106); therefore, even though a connection between them 
was evidenced, it was relatively weak. In summary, we found that years of 1961 and 1974 were 
not only the main trend change points of the dependent variable (water volume change), as seen in 
Figure 3, but also the years with the highest repetition rate and probability among the other 
variables, ranked the second and first, respectively. We used the stage division as the basis for 
further trend and statistical analyses. 


3.1.2 Trend analysis and statistical characteristics by periods 


We conducted trend (Table 1) and statistical (Table 2) analyses for each variable over three 
periods (1950-1960, 1961-1973, and 1974-2016). In periods one and two (1950-1960 and 
1961-1973), all variables were stationary, but most showed distinctive statistical characteristics 
(including mean, minimum, and maximum values). Among them, the values for the water volume 
change, K-S variable, and K-K variable became larger during the transition from period one to 
period two, while variables Syr Down_R and Samanbay_R displayed the opposite behaviour. In 
period three (1974-2016), the Syr Down_R, K-K, and Kelif_R variables remained flat, but their 
specific values became smaller, larger, and smaller, respectively, compared to those in period two. 
The water volume change showed a significant downward trend with a P-value<0.01, a slope 
value of —1.032, and a mean value of 17.719 km?/a. The Samanbay_R variable showed a weak 
downward trend with a P-value<0.05, a slope value of —0.175, and a mean value of 7.469 km’. 
The K-S variable showed a weak upward trend with a P-value<0.05, a slope value of 0.175, and a 
mean value of 51.076 km’. In particular, the K-K variable was further discussed in turn for other 
three periods: 1950-1958, 1959-1990, and 1991-2016, because of its distinctive trend 
characteristics. Specifically, it started with a non-stationary sequence in period one (1950—1958), 
showing a significant downward trend with a slope value of —1.746 and a mean value of 0.026 km’. 
Then, it transitioned to another non-stationary sequence in period two (1959-1990), showing a 
significant upward trend with a slope value of 0.341 and a mean value of 11.709 km’. Finally, it 
transitioned to a stationary sequence in period three (1974-2016), with a mean value of 16.973 


km’, 


Table 1 Trend analysis of the variables by different periods and their slope values 


1950-1960 1961-1973 1974-2016 
Variable 
Trend Slope P Trend Slope Ẹ Trend Slope P 
WVC > - - > - - l —1.032 ss 
Syr_Down_ R — - - — - - — - - 
Samanbay_R <> - - > - - Li —0.175 we 
K-S => - - — - - tT 0.175 ae 
K-K =# = = = m = = x m 
Kelif R > - - > - - > - - 
1950-1958 1959-1990 1991-2016 
Variable 
Trend Slope P Trend Slope P Trend Slope P 
K-K l —1.746 me tT 0.341 ie > - - 


Note: —, horizontal cyan arrow represents no significant upward or downward trend; Î, upward red arrow represents an upward trend; 
|, downward purple arrow represents a downward trend; WVC represents the water volume change in the Aral Sea; Syr_Down_R, 
runoff from the downstream of the Syr Darya River; Samanbay_R, runoff from the downstream of the Amu Darya River; K-S, water 
consumption in the middle-lower reaches of the Amu Darya River; K-K, water consumption in the middle-upper reaches of the Amu 
Darya River; Kelif_R, upstream-originating flow of the Amu Darya River. ** indicates significant level at P<0.05 level; *** indicates 
significant level at P<0.01 level. -, no data. 
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Table 2 Statistical characteristics of the variables by different periods 


1950-1960 1961-1973 1974-2016 
Variable 
Mean Min Max Mean Min Max Mean Min Max 
WVC (km*/a) -3.103 -14.076 6.724 15.212 —6.112 30.570 17.719 —6.033 48.968 


Syr_Down_R (km) 16.127 9.500 21.100 6.400 3.200 10.600 0.908 0.200 2.100 
Samanbay R(km*) 45.994 31.005 55.400 35.059 21.818 71.067 7.469 0.117 24.196 


K-S (km*) 12.550 3.145 27.539 23.486 -12.523 36.727 51.076 34.348 58.427 
K-K (km°) 0.630 8.200 9.756 7.948 3.527 14.836 16.295 10.224 26.723 
Kelif_R (km) 67.182 49.600 80.000 63.485 50.000 96.300 59.837 35.300 83.500 
1950-1958 1959-1990 1991-2016 
Variable 
Mean Min Max Mean Min Max Mean Min Max 
K-K (km’) 0.026 8.200 9.756 11.709 2.375 22.005 16.973 11.750 26.723 


Note: Min indicates minimum; Max indicates maximum. 


3.2 Relationship between the water volume change in the Aral Sea and its recharge runoff 


3.2.1 Non-linear cross-correlation based on the MF-DCCA 


In this study, the MF-DCCA method was used to test the long-range cross-correlation of the water 
volume change with the variables of Samanbay_R, K-S, K-K, Kelif_R, Kerki_R, and water area, 
by obtaining the Hurst index at different g-values. As shown in Figure 4, where the horizontal axis 
is the q-value and the vertical axis is the Hurst index value, there was a large difference in the 
Hurst index obtained by the water volume change and each variable at different g-values, 
indicating multifractal characteristics in their mutual correlation. When g-values successively 
increased from —5 to 5, the Hurst index was different in the non-linear functions, indicating that 
there were different power-law relations between the water volume change and other variables. By 
comparing the variation range of these cross-correlation indices, we found that the 
cross-correlations between the water volume change and water area had the highest multifractal 
degree, while those between the water volume change and Kelif R variable had the lowest 
multifractal degree. This indicates that the correlation between the water volume change and water 
area was the most complex and had the strongest linkage, while the correlation between the water 
volume change and Kelif_R variable was the simplest and weakest. This is attributed to the 
relationship between the water volume change and water area being both non-linear, as determined 
by topography, and also representing a joint response to the net water outflow from evaporation. 

The lowest Hurst index values between the water volume change and Samanbay_R variable, 
between the water volume change and K-S variable, and between the water volume change and 
K-K variable were 0.66, 0.66, and 0.64, respectively. These values were all greater than 0.50, 
suggesting that there exists a long-range cross-correlation between them, regardless of the 
fluctuation magnitude; that is, when a sequence had a certain trend, another sequence was likely 
to show a synergistic trend. Of these, the variable pairs of the water volume change-Samanbay_R 
and water volume change—K-S not only had the smallest Hurst index values greater than 0.50, but 
also showed the same overlapping curve in Figure 4, proving that the Samanbay_R and K-S 
variables have mostly consistent synergistic effects on the water volume change. 

However, the cross-correlation of the water volume change with the variables of Kelif_R, water 
area, and Kerki R denoted a power-law relation that was completely different from that of the 
previous variables. When g-value<0, each cross-correlation index was smaller than 0.50, showing 
that the cross-correlation at smaller fluctuations was anti-persistent; that is, smaller changes 
between the water volume change and other variables were more likely to change in the opposite 
direction. When the Kelif_R variable showed a big upward trend, for example, there was a very 
likely sharp downward tendency in the water volume change, which was also consistent with the 
large increase in the upstream water supply and the slowing process of the Aral Sea shrinkage. 
When q-value>0, each cross-correlation index was larger than 0.50, showing that the 
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cross-correlation at larger fluctuations had a relatively long-range persistence. Namely, the large 
fluctuations in the water volume change would be affected by the long-term effects of the 
Kelif R, Kerki_R, and water area, which were more likely to change in the same direction. For 
example, when water area showed a small upward (or downward) trend, the water volume change 
was also likely to demonstrate such a trend. 

Additionally, we found that the water volume change and water area had the strongest and most 
complex long-range correlation. The long-range correlation of the water volume change with the 
Kelif R variable, the Kerki R variable, and the Samanbay_R variable, also successively 
strengthened. This may have been related to the growing geographical proximity to the Aral Sea. 
The persistent influence on the water volume change by the K-K variable and K-S variable showed 
the opposite state under different fluctuation amplitudes, which may related to the mechanism of 
their influence on the Aral Sea shrinkage. The water consumption of different reaches could not 
only directly affect the replenishment of water into the Aral Sea, but also indirectly impact the 
replenishment of water into the Aral Sea through the influence of human activities (e.g., irrigation). 
Consequently, various mechanisms produced different long-range effects. 
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Fig. 4 Long-range cross-correlations of the water volume change (WVC) with different variables. q-value is the 
order of the fluctuation function that can be changed to examine different characteristics of the data. 


3.2.2 Quantitative analysis 


The model simulation results are presented in Tables 3 and 4. For Model I, by constructing the 
relationships of the water volume change with water area, Samanbay_R, and Syr_Down_R, and 
combining them with the theoretical capacity of the Aral Sea, we obtained quantitative 
relationships between the water volume of the Aral Sea and its recharge runoff during 1950-1985. 
The model was statistically significant, and each coefficient was reliable. As shown by the 
standardised coefficients, the relative importance of the net evaporation (represented by water 
area) and the recharge runoff of the Aral Sea (represented by the Samanbay_R and Syr_Down_R) 
to the Aral Sea shrinkage was 0.170, 0.741, and 0.401, respectively. Model II obtained 
quantitative relations between the water volume and the runoff in the Amu Darya River during 
1986-2016 by constructing the relations of the water volume change with the water area and 
Samanbay_R. The model was statistically significant, and each coefficient was reliable. As can be 
seen from the standardised coefficients, the relative importance of the water area (net 
evaporation) and Samanbay_ R to the Aral Sea shrinkage amounted to 0.647 and 0.690, 
respectively. Model III showed a quantitative relationship of the water volume with the runoff in 
the Amu Darya River and water consumption of the Amu Darya River during the 1986-2016 
period by constructing the relationship of the water volume change with the water area, K-K, K-S, 
Kelif R and THC_V. Overall, the model fitted well but the coefficients for the K-K were not 
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statistically significant, unlike the other coefficients. Therefore, the direct effect of the K-K 
during the 1986-2016 period was insignificant. We removed the K-K variable, continued the 
construction of the quantitative models, and obtained Model IV. Model IV was overall statistically 
significant, and its respective coefficients were reliable. The standardised coefficients showed that 
the relative importance of the net evaporation, water consumption of the Amu Darya River 
(excluding the reservoir water), THC_V, and Kelif_R to the shrinkage of the Aral Sea were 0.934, 
0.705, 0.446, and 0.054, respectively. 


Table 3 Specific structure and parameter estimates of Model I and Model II 


Model I (1950-1985) Model II (1986-2016) 
WVC=a,x WA-f, xR, -£,xR, WVC =a, x WA- pa xR, 
2.68e —9 x WV" + 3.28 x WV? +1.59e -3 x 
WV? -0.41x WV +11.54. 
0.08 x WV — 20.86, WV €(961,1006] x + 
Formula : WV €(79,422] 
WA= 4 —1.68e - 05x WV? + 0.05 x WV + WA = RE E E E E T, 
22.52, WV €(422,961] ee ee ee ee 
WV —3.75 x WV + 51.31, 
WV € [43,79] 
GOF R°=0.9212 Adjusted R?=0.9164 R°=0.9326 Adjusted R?=0.9303 
Model I: parameter estimate Model II: parameter estimate 
ay Bu p2 a2 bn 
Coef. 0.837" 0.772" -1.217™" 0.864" -0.964"" 
Std coef. 0.170" 0.741" —0.401™ 0.647" —0.690"™ 


Note: GOF represents the goodness-of-fit of the model; Coef. represents the unstandardized coefficients for the multiple regression 
equations in the models; and Std coef. represents the standardised coefficients for the multiple regression equations in the models. WA 
represents the water area in the Aral Sea; R; represents the recharge runoff from the Amu Darya River; R, denotes the recharge runoff 


from the Syr Darya River; and a, a2, J11, Bi2, and £2 are the coefficients. ** indicates significant at P<0.05 level, and *** indicates 
significant at P<0.01 level. 
Table 4 Specific structure and parameter estimates of Model III and Model IV 
Model III (1986-2016) Model IV (1986-2016) 
WVC=a,x WA-7xQ-7 x0, =7,xU -yx RV WVC=a,x WA- xC, -7 xU -y,x RV 
2.68e —9 x WV" + 3.28 x WV? +1.59e -3x 2.68e —9 x WV" +3.28x WV? +1.59e -3x 
WV? - 0.41 x WV 411.54, WV? - 0.41x WV 411.54, 
Formula WV €(79, 422] WV €(79, 422] 
WA= WA= 
3.97e — 6x WV" —1.04e —3 x WV? + 0.10 x 3.97e — 6 x WV" —1.04e -3 x WV? + 0.10 x 
WV? -3.75x WV + 51.31, WV? -3.75 x WV + 51.31, 
WV €[43,79] WV € [43,79] 
GOF R°=0.9074 Adjusted R7=0.8931 R°=0.9041 Adjusted R?=0.8934 
Model III: parameter estimate Model IV: parameter estimate 
a3 yı y2 y3 Y4 a3 y2 y3 Ya 
Coef. 0.765** 0.201 0.360"  —-0.348"" 0.204" 0.768% 0.382% 0.323" 0.246°* 
Std coef. 0.934" 0.010 0.705"  —0.054 0.446" 0.936 0.713" —0.045 0.454 


Note: Cı represents the water consumption in the middle-upper reaches of the Amu Darya River (K-K); C, represents the water 
consumption in the middle-lower reaches of the Amu Darya River (K-S); U represents the upstream-originating flow of the Amu Darya 
River (Kelif_R); RV represents the water storage in the reservoir of the Tuyamuyun Hydro Complex (THC_V); and a3, yı, y2, y3, and y4 
are the coefficients. 


In summary, during the 1950-1985 period, the Samanbay_R was the dominant factor affecting 
the Aral Sea shrinkage; the second was the Syr Down_R, followed by the net evaporation. 
During 1986-2016, the water loss by net evaporation and the K-S were almost equally important 
to the water volume, even though net evaporation was far more influential than water use along 
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the river sections of the Amu Darya River. The importance of the variables affecting the Aral Sea 
followed an order of K-S>THC_V>Kelif_R>K-K. 


4 Discussion 


4.1 Historical evolution of the Aral Sea shrinkage and its causes 


Based on the above analysis of trend change points, we discussed the evolutionary process of the 
Aral Sea in different time periods, as well as the changes in other hydrological process associated 
with it. By combining with key background events in the Aral Sea region during 1950-2016 
period (Table 5), we explored the specific causes of the Aral Sea shrinkage. 

From 1950 to 1960, the water volume change tended to be flat and showed a small negative 
value, implying a weak increasing trend in the water volume of the Aral Sea. However, this state 
was interrupted in 1961, after which the Aral Sea started to shrink, with a shrinkage rate 
fluctuating around 15.212 km?/a. Correspondingly, the Samanbay_R and Syr Down_R underwent 
abrupt changes in 1961, with a sudden decrease in runoff. The former had an average annual 
runoff value of 35.059 km? during 1961-1973, showing a decline of 10.935 km? from the 
previous time period, whereas the latter had an average annual value of 6.400 km’, representing a 
decrease of 9.727 km? from the former time period. We believe that the abrupt decrease in the 
downstream of the recharging rivers during the 1960s caused the Aral Sea to shrink because the 
received recharge amount seemed insufficient to compensate the evaporation. This result is 
attributed to the large-scale development of irrigated agriculture in the Aral Sea region starting in 
the late 1950s (Aldaya et al., 2010; Bekchanov et al., 2012; Tischbein et al., 2012), when the 
Soviet government ordered a specialisation in cotton production. Consequently, the construction 
of the gigantic Karakum Canal began in 1954 and started working in 1956, drawing water from 
the Amu Darya River to Turkmenistan (Micklin and Aladinet al., 2008). After 1960, the amount 
of water used for irrigation increased (Spoor, 1993; Abduraupov et al., 2022), as did the irrigated 
area from approximately 4.510*104 km? in 1960 to about 5.150*10* km? in 1970. 

The year 1974 represented the most significant change point in change point detections for all 
time series. During this year, the Aral Sea shrank at a rate of up to 48.968 km*/a, which was the 
highest value measured in this study. Simultaneously, the Samanbay R and Syr Down R 
plummeted to 9.049 and 1.300 km’, respectively. The continued significant decrease in the flow 
of the two recharging rivers is believed to have caused a sudden and accelerated shrinkage of the 
Aral Sea, with the Samanbay_R and Syr Down R dropping by up to 26.010 and 5.100 km’, 
respectively, compared to the average values registered during the 1961-1973 period. During the 
1970s, irrigation technologies were greatly developed, and a large number of irrigation facilities 
were built in the Aral Sea region (Abdullaev, 2004). For example, the Karshi Canal and 
Amu-Bukhara Canal, two of the most important water diversions along the Amu Darya River, 
were put into operation in 1973 and 1974, respectively (corresponding to the change point year of 
1973), with respective maximum flows of 160.000 and 360.000 m?/s and average annual 
diversions of 4.000 and 4.500 km? (Schliiter et al., 2005). This suggests that there is a qualitative 
improvement in the efficiency of the local use of water resources. The large amount of water used 
for agriculture in the 1970s together with new land development schemes in the Hunger, Karshi, 
Dzhizak, and Sherabad steppes, the Karakum Canal, and Kyzylkum Canal (Dukhovny et al., 2007) 
led to the above-mentioned phenomenon, that is, a sharp decrease in the Samanbay R and 
Syr_Down_R and a high rate of the shrinkage of the Aral Sea. 

Following, the Syr Down R remained steady during the 1974-1985 period, while the 
Samanbay_R showed a weak and stable decreasing trend (a slope value of —0.175 km?/a during 
the period 1974-2016), with a very small average value that can be explained by human activities 
at that time. The irrigated area in the Aral Sea region continued to grow in the 1980s but stagnated 
in the late 1990s, after which it stabilized at 8.000x10* km* (FAO, 2020). This is due to the 
collapse of the former Soviet Union in 1991, when the construction of a large number of irrigation 
facilities ceased, followed by a reduction in the cotton area and a decline in cotton production 
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(Spoor, 1998). Although the irrigated area ceased to grow significantly after the 1990s, the 
cropping structure began to shift, with wheat gradually replacing cotton as the main local crop 
(Abdullaev et al., 2005; Kariyeva and van Leeuwen, 2012; de Beurs et al., 2015), so agricultural 
water use continued to grow. For example, in Uzbekistan, the share of cotton in the total irrigated 
area significantly declined during the 1990s, while the share of grain crops increased. However, 
the shrinkage rate of the Aral Sea has been continuously decreasing, with a slope value of —1.032 
km/a, implying that the Aral Sea continued to shrink at a gradually slower rate. This is because 
during this period, as the shrinking area of the Aral Sea was decreasing, evaporation also 
gradually declined, leading to a gradual drop in the shrinkage rate and a continuous advance 
towards the Aral Sea, achieving re-equilibrium (i.e., the evaporation and recharge were nearly 
equal). At this stage, the importance of net evaporation for the evolution of the Aral Sea began to 
emerge, though it has been ignored in most existing studies (Yang et al., 2020b; Li et al., 2021; Su 
et al., 2021; Shi et al., 2022). 

Because the Samanbay_R was determined by the amount of incoming upstream-originating 
flow and the amount of water consumed by each river section, the Kelif_R, K-K, and K-S 
variables were considered in this study. Considering the described results, the trend of the Kelif R 
was relatively stable, without evident upward or downward trends. The K-S variable kept rising 
during three periods (1950-1961, 1962-1973, and 1974-2016), with a rising proportion in the 
total water consumption. Meanwhile, the magnitude of the K-K variable was much smaller than 
that of the K-S variable, and the change point years of the K-K variable were completely different 
from those of the water volume change variable, meaning that its influence on the Aral Sea was 
not as strong as that of the K-S variable. Therefore, we believe that the K-S was the main reason 
for the variation of the Samanbay_R. Moreover, the results of the change point and long-range 
cross-correlation analyses were consistent with the Samanbay_R variable, with the same change 
points (in 1961 and 1974, respectively) and long-range correlation curve. 


Table 5 List of the background events in the Area Sea region during 1950-2016 period 
1950-1960 1961-1973 1974-1991 1992-2016 
(1) The Soviet Union 


C1) The Soviet government (1) The Karshi Canal began 


ordered the large-scale (1) From 1960 to 1970, its operation in 1973 dissolved in 1991. 
cotton production in the the irrigated area in the (2) Te A mü- Bukhara Canal (2) Around 1992, 
late 1950s. Aral Sea region significant changes 


Events was put into use in 1974. 


(2) The construction of the experienced an (3) The cott F occurred in the cropping 
Karakum Canal began in approximate increase; i PEREO WNE Area structure, with the cotton 
: 103 km2 in the Aral Sea region reached x 
1954, and was put into of 6.400* 10° km’. : area decreasing and the area 
we a stable state in the 1980s. : : : 
operation in 1956. of grain crops increasing. 


4.2 Exploring the response of the Aral Sea to human activities by scenario analysis 


Using various historical water use scenarios, the evolution of the Aral Sea under different 
scenarios was simulated to better understand the impact of different water use policies. According 
to the water use analysis and purpose of this study, combined with existing studies (Froebrich and 
Kayumov, 2004; Sun et al., 2019; Yuldashev et al., 2020), it should be possible to increase the 
water coming from the lower reaches of the Amu Darya River to different degrees of the water 
level in 1961-1973 by increasing the efficiency of agricultural water use and reducing it in the 
middle and lower reaches, thereby reducing ineffective evaporation from reservoirs and wetlands. 
Thus, we set up three scenarios: Scenario 1, after 1974, with the Samanbay_R without a cliff drop 
and maintaining the flow level of 1961—1973 in the previous period (i.e., the water consumption 
of the river section (upstream of the Samanbay hydrological station) maintained the flow level of 
1961-1973 during 1974-2016); Scenario 2, compared to the flow level during the period 
1974-2016 under Scenario 1 with a 100% recovery, it is assumed that the Samanbay_R recovered 
by 50% after 1974 (1974-2016); and Scenario 3, assuming that no THC reservoir (in-channel) 
was historically built (it was actually constructed in 1983), namely, the water consumption of the 
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river section (upstream of the Samanbay hydrological station) continued to maintain the same 
level during the historical period of 1974-1982. Scenarios 1 and 2 were assumptions about 
whether agriculture expanded along the Amu Darya River after 1974 and to what extent. 
Contrastingly, Scenario 3 was designed to explore whether the reservoir water storage had an 
effect on the Aral Sea. The Samanbay_R during 1974-2016 period was 100%, 50%, and 
approximately 30% of the flow level during 1961—1973 period under the Scenario 1, Scenario 2, 
and Scenario 3, respectively (Fig. 5). 
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Fig. 5 Comparison between historical runoff values in the downstream of the Amu Darya River and simulated 
runoff values under different scenarios 


Through the joint simulation of Model I and Model II shown above, we obtained the 
evolutionary process of the water volume in the Aral Sea under the three scenarios (Fig. 6). 
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Fig. 6 Comparison between historical water volume values in the Aral Sea and simulated water volume values 
under different scenarios 


Under Scenario 1, after 1974, the Samanbay_R remained consistent with the previous level 
(average annual value of 35.060 km? during 1961-1973), while the total annual water 
consumption of the Amu Darya River was approximately 24.778 km’. In this case, during the 
evolution of the Aral Sea, demonstrated by the blue line, the water volume decreased slowly, from 
843.023 km? in 1973 to 666.693 km? in 1985 (221.987 km? higher than the historical value of 
444.706 km), until it fell to 462.185 km? in 2016, well above the historical value (54.561 km’). 
There was no split since then. Under Scenario 2, after 1974, when the Samanbay_R was 
maintained at an average annual level of approximately 21.264 km, the total annual water 
consumption of the Amu Darya River was approximately 38.573 km’. In this case, the water 
volume of the Aral Sea gradually decreased from 843.023 km? in 1973 to 383.803 km? in 1996, 
leading the Aral Sea dividing into the Big Aral Sea and Small Aral Sea, with a 10-year delay 
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compared with history (1986). Additionally, by 2016, water volume decreased to 183.552 km’. 
Under Scenario 3, the THC reservoir was no longer constructed in 1983 and the Samanbay_R 
maintained the average annual level of approximately 10.351 km? prior to 1983; that is, the total 
annual water consumption of the Amu Darya River measured approximately 46.916 km’. In this 
case, the water volume rapidly decreased from 843.023 km? in 1973 to 394.062 km? in 1986, 
when the Aral Sea was divided into the Big Aral Sea and Small Aral Sea. This time point was 
consistent with history, showing that the suspension of the THC reservoir in 1983 could not 
change the evolution of the Aral Sea, as indicated in Figure 6. After the division of the Aral Sea in 
1986, from 1986 to 2016, the water volume of the Aral Sea was slightly larger than the historical 
water amount (approximately 33.583 km*), and the water volume fell to 83.251 km? until in 2016. 


5 Conclusions 


The famous Aral Sea crisis in the arid zone of Central Asia, which poses huge pressures on the 
local ecology and sustainable development, is mainly attributed to the excessive exploitation of 
water resources in the Amu Darya River by human activities. In this study, we applied the BEAST 
model to detect the trend change points of involved variables, discussed the non-linear 
relationships between the Aral Sea's water volume and Amu Darya River's runoff using the 
MF-DCCA method, and provided their quantitative correlations. On this basis, we explored the 
evolution of the Aral Sea during the 1950-2016 period and analysed the causes of its shrinkage 
over different time periods. Additionally, we explored how the Amu Darya River affected the Aral 
Sea both historically and under different scenarios. The main conclusions are as follows: 

(1) During the 1950-2016 period, the water volume change of the Aral Sea had two significant 
trend change points in 1961 and 1974. Simultaneously, other hydrological variables linked to the 
Aral Sea, of which, the Samanbay_R, Syr Down_R, and K-S, showed the same trend change 
points in 1961 and 1974, the Kelif R variable had only one weak trend change point in 1961, 
whereas the K-K variable had two completely different trend change points, in 1959 and 1991. 
These change points, corresponding to significant changes in time series data, help reveal the state 
of the Aral Sea and the identified variables during different periods. The years of 1961 and 1974, 
frequently appeared in various related variables, implying their strong correlation between each 
other. 

(2) The long-range cross-correlation between the water volume change of the Aral Sea and the 
related variables (Samanbay_R, Kerki_R, Kelif R, K-S, K-K, and water area) showed 
multifractal characteristics. The water volume change showed long-range cross-correlations with 
the Samanbay_R, K-S, and K-K variables consistently, regardless of the fluctuation magnitude. 
The water volume change only had long-range cross-correlations at larger fluctuations with the 
Kerki_R, Kelif R, and water area variables. Finally, the Samanbay_R and K-S variables had 
mostly consistent synergistic effects on the water volume change. 

(3) Before 1974, the Samanbay_R was the dominant factor affecting the Aral Sea shrinkage, 
followed by the Syr Down_R and the net evaporation. From 1974, the importance of net 
evaporation began to appear, and after 1986, its impact on the Aral Sea was almost as important 
as the recharge water of Samanbay_R, even though it was far more influential than water use 
along the river sections of the Amu Darya River. In summary, the decrease in the Samanbay_R 
was the main cause of the Aral Sea shrinkage, and the K-S was the main reason for the variation 
in the Samanbay_R. The order of importance of the variables in the Amu Darya River system 
affecting the Aral Sea shrinkage is as follows: K-S>THC_V>Kelif_R>K-K. 

(4) This study evaluates the evolution of the Aral Sea under different scenarios and compares it 
with the historical evolution process. It was found that when the Samanbay_R during the 
1974-2016 period recovers by 50% of the 1961-1973 level, the Aral Sea would split in 1996, 
which is a decade later than the historical split in 1986. Thus, it is possible to slow down the 
retreat of the Aral Sea and stabilize its water area and water volume in an ecologically restored 
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state, when recovering the Samanbay_R after 1974 to the 1961-1973 level, by increasing the 
efficiency of agricultural water use, reducing agricultural water use in the middle and lower 
reaches, and reducing ineffective evaporation from reservoirs and wetlands. 

This study examined the relationship between the Aral Sea shrinkage and the runoff of the Amu 
Darya River in detail and offers new perspectives on the causes of shrinkage in the Aral Sea. 
However, further improvement is required, and more human activity factors and human policies 
on water and land use should be taken into consideration in the future. 
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